import matplotlib.pyplot as plt
from mpmath import  mp,mpf,exp,pi,sqrt

# 显示3个不同版本的1/sqrt(x)的 Remez 近似多项式的误差
# 设置多精度浮点数的精度
mp.dps = 50  # 可以根据需要调整精度

# 定义一个使用mpf的函数
def my_invsqrt_error_v1(x):
    #C0 = 1.1033542890963093621648566115027039
    #C1 = -0.1666666666666666666666666666666667
    C0 = mpf("1.1033542890963093621648566115027039")
    C1 = mpf("-0.1666666666666666666666666666666667")
    
    r = C0 + C1*x
    return r- (1.0/sqrt(x))

def my_invsqrt_error_v2(x):
# C0 = 1.3630036810814256614681164551181137
# C1 = -0.4372649921698052571543343652873008
# C2 = 0.0564892705109514684020846421294280
    C0 = mpf("1.3630036810814256614681164551181137")
    C1 = mpf("-0.4372649921698052571543343652873008")
    C2 = mpf("0.0564892705109514684020846421294280")
    
    r = C0 + C1*x  + C2*(x**2)
    return r- (1.0/sqrt(x))

def my_invsqrt_error_v3(x):
# C0 = 1.5800398739333915918140553579337426
# C1 = -0.7743678147619097758162146830650696
# C2 = 0.2102255238402672409306734208569276
# C3 = -0.0211155462431472902620866232326779
    C0 = mpf("1.5800398739333915918140553579337426")
    C1 = mpf("-0.7743678147619097758162146830650696")
    C2 = mpf("0.2102255238402672409306734208569276")
    C3 = mpf("-0.0211155462431472902620866232326779")
    
    r = C0 + C1*x  + C2*(x**2) + C3*(x**3)
    return r- (1.0/sqrt(x))

def plot_exp_error(mode):

    k=(1.0/100)
    x_values = [ mpf(i)*k for i in range(100,401)]  
    
    # 计算对应的y值
    if mode==1:
        y_values = [my_invsqrt_error_v1(x) for x in x_values]
    elif mode==2:
        y_values = [my_invsqrt_error_v2(x) for x in x_values]
    else:
        y_values = [my_invsqrt_error_v3(x) for x in x_values]
   
    # 绘制图像
    plt.plot(x_values, y_values)
    plt.xlabel('x')
    plt.ylabel('y')
    
    if mode==1:
        plt.title('r(x)-invsqrt(x)')
    elif mode==2:
        plt.title('r2(x)-invsqrt(x)')
    else:
        plt.title('r3(x)-invsqrt(x)')
    

    plt.grid(True)
    plt.show()


plot_exp_error(1)
plot_exp_error(2)
plot_exp_error(3)
